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Numerous  experimental  investigations  clearly  established  that  when  soda-lime  glass  is  subjected  to  suf¬ 
ficiently  high  axial-stress/pressure,  it  displays  a  nonlinear  mechanical  response  and  deformation  irrevers¬ 
ibility  (inelasticity).  This  portion  of  the  material  behavior  is  often  neglected  in  material  models  for  glass 
which  tend  to  focus  on  the  damage  and  fracture  phenomena  of  the  material.  However,  material  nonlin¬ 
earity/inelasticity  can,  in  principle,  have  a  profound  effect  on  wave/shock  propagation  phenomena  and 
processes  (e.g.  spall  fracture). 

Within  the  present  work,  the  effect  of  material  nonlinearity  and  inelastic  behavior  on  the  dynamic 
response  (including  spallation)  of  soda-lime  glass  is  studied  under  symmetric  flyer-plate  loading  condi¬ 
tions  using  computational  methods  and  tools.  Material  nonlinearity  and  deformation  irreversibility  are 
modeled  in  two  different  ways:  (a)  as  a  non-linear  elastic  material  response  with  no  deformation  irre¬ 
versibility;  and  (b)  as  a  linear-elastic,  volumetrically-plastic  deformation  response.  Incorporation  of  non¬ 
linearity  and  inelasticity  phenomena  into  a  continuum-level  material  model  for  soda-lime  glass  recently 
developed  by  the  authors  revealed  that  while  these  phenomena  do  not  measurably  affect  spall  resistance 
(as  measured  by  a  minimum  flyer-plate  velocity  resulting  in  spallation),  they  provide  beneficial  linear- 
momentum/kinetic  energy  reduction  effects. 

©  201 1  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Contemporary  transparent  armor  is  typically  composed  of  a 
system  of  materials  designed  to  meet  a  set  of  optical  transparency 
constraints  while  providing  a  compulsory  level  of  protection 
against  blast  and  ballistic/fragment  impacts.  This  class  of  protective 
materials  is  employed  in  varied  applications  from  personal 
protective  visors  for  combat/non-combat  usage  (e.g.  riot  control 
or  explosive  ordinance  disposal)  to  vehicle/structure  transparent- 
armor  systems  (vehicle  windows  designed  to  protect  occupants 
and  on-board  instruments/sensors  from  projectiles  and/or 
fragment  impacts  during  combat/terrorist  attacks  or  other  hostile 
conflicts).  Recent  engagements  of  the  US  military  forces  in  the 
Operation  Iraqi  Freedom  and  Operation  Enduring  Freedom 
(Afghanistan)  have  highlighted  the  critical  importance  of  transpar¬ 
ent  armor.  Persistent  escalations  in  the  number  and  variety  of 
threats  along  with  constricting  military  budget  have  greatly 
increased  the  need  for  rapidly-deployable,  threat-specific, 
weight-/cost-performance-optimized  transparent  armor  and  armor 


*  Corresponding  author.  Address:  241  Engineering  Innovation  Building,  Clemson 
University,  Clemson,  SC  29634-0921,  United  States.  Tel.:  +1  864  656  5639;  fax:  +1 
864  656  4435. 

E-mail  address:  gmica@clemson.edu  (M.  Grujicic). 

0261-3069/$  -  see  front  matter  ©  2011  Elsevier  Ltd.  All  rights  reserved. 
doi:10.1016/j.matdes.201 1.08.031 


systems.  These  threats  to  US  military’s  safety  have  motivated  great 
amount  of  research  effort  within  US  and  abroad  focused  on  acceler¬ 
ating  the  development  of  novel/enhanced  transparent  armor 
systems. 

Historically,  transparent  armors  have  been  based  on  monolithic 
glass  or,  more  recently,  transparent-elastomer  inter-layered  glass 
laminates.  Numerous  advancements  within  the  class  of  transpar¬ 
ent-armor  materials  and  novel  technologies  have  been  reported 
within  open  literature,  of  which  the  following  have  received  the 
most  attention:  crystalline  ceramics  (e.g.  aluminum-oxinitride  spi¬ 
nel,  A10N  [  1  ]),  novel  transparent  polymer  materials  (e.g.  transparent 
nylon  [2]),  new  interlayer/adhesive  materials  (e.g.  polyurethane 
bonding  layers  [1]),  and  new  laminate  designs  e.g.  [3].  In  the  face 
of  increasing  demands  for  improvements  in  ballistic-protection  per¬ 
formance  of  transparent  armor,  and  the  concomitant  requirements 
for  improved  performance  to  weight  ratios  that  call  for  the  use  of 
new  transparent  materials  (e.g.  transparent  crystalline  ceramics), 
advanced  transparent  polymeric  materials  and  advanced  technolo¬ 
gies  (e.g.  multi-material  functionally-graded  laminated  transparent 
armor),  glass  (as  well  as  glass  ceramics)  continues  to  retain  its  role  as 
a  popular  material  choice  in  ground-vehicle  transparent  armor 
applications.  The  main  reason  for  the  continued  use  of  glass  in  trans- 
parent-armor  applications  have  been  discussed  in  our  recent  work 
[4,5]. 
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It  is  well  understood  that  glass  exhibits  markedly  different 
behavior  under  quasi-static  (i.e.  low  deformation-rate)  and 
dynamic  (i.e.  high  deformation-rate)  loading  conditions  e.g. 
[6-12].  Main  differences  in  glass  behavior  within  the  two  loading 
rate  ranges  have  been  discussed  in  our  recent  work  [4,5].  Further¬ 
more,  the  physics  of  glass  behavior  in  these  two  regimes  has  been 
utilized  in  Refs.  [4,5]  to  construct  a  continuum-type  material 
model  for  glass. 

A  comprehensive  literature  review  carried  out  as  part  of  the 
present  work  revealed  that  the  mechanical  behavior  of  glass  is 
modeled  predominantly  using  three  distinct  approaches:  (a) 
molecular-modeling  methods  [13-19];  (b)  continuum-level  finite 
element  material  approximations  [20-24,26,27],  and  (c)  finite- 
element  modeling  of  explicit  cracks  [25].  A  brief  overview  and 
the  main  findings  for  each  of  these  three  approaches  can  be  found 
in  our  recent  work  [4,5].  The  present  work  deals  with  the  contin¬ 
uum-level  glass  models  which  are  generally  used  in  large  scale 
simulations  of  the  ballistic  performance  of  transparent  armor 
systems. 

As  mentioned  earlier,  a  continuum-level  material  model  for 
glass  has  been  proposed  in  our  recent  work  [4,5].  The  main  objec¬ 
tive  of  the  present  work  is  to  extend  this  model  in  order  to  include 
the  phenomena  of  material  non-linearity  and  deformation- 
irreversibility.  To  assess  the  contribution  of  these  phenomena  to 
the  dynamic  behavior  of  glass,  computational  flyer-plate  experi¬ 
ments  are  carried  out  and  the  spallation  fracture  process 
monitored. 

The  organization  of  the  paper  is  as  follows:  An  overview  of  the 
soda-lime  glass  material  model  of  Grujicic  et  al.  [4,5]  is  provided  in 
Section  2.  The  incorporation  of  the  material-nonlinearity  and 
deformation  inelasticity  effects  into  the  material  model  of  Grujicic 
et  al.  [4,5]  is  presented  in  Section  3.  Details  of  a  transient  non¬ 
linear  dynamics  computational  analysis  of  a  flyer-plate  impact 
experiment  used  to  validate  the  enhancements  to  the  material 
model  for  soda-lime  ballistic  glass  are  discussed  in  Section  4. 
The  main  results  obtained  in  the  present  work  are  presented  and 
discussed  in  Section  5.  The  key  conclusions  resulted  from  the 
present  work  are  summarized  in  Section  6. 


2.  Soda-lime  glass  model  of  Grujicic  et  al.  [4,5] 

This  section  focuses  on:  (a)  the  main  features  of  the  soda-lime 
glass  material  model  developed  by  Grujicic  et  al.  [4,5];  (b)  the  pro¬ 
cedure  used  to  derive  the  fundamental  equations  which  govern  the 
material  mechanical  response;  and  (c)  the  approach  used  to  gather 
the  necessary  experimental  data  and  carry-out  material-model 
parameterization. 

2.1.  Damage/failure  physics  of  soda-lime  glass 

The  soda-lime  glass  model  of  Grujicic  et  al.  [4,5],  tries,  in  a  com- 
putationally-efficient  manner,  to  incorporate  the  experimentally 
well-documented  differences  in  the  behavior  of  glass  when 
subjected  to  low  deformation-rate  (quasi-static)  and  high  defor¬ 
mation-rate  (dynamic)  loading  conditions.  Specifically,  it  is  well- 
established  that,  under  quasi-static  loading  conditions,  failure  of 
glass  is  associated  with  the  formation/propagation  of  few  cracks 
and  the  complete  fracture  yields  few  fragments.  On  the  other  hand, 
under  dynamic  loading  conditions,  glass  failure/fracture  is  pre¬ 
ceded  by  extensive  damage  (associated  with  the  formation  of 
numerous  micron  and  submicron-size  cracks).  The  complete  frac¬ 
ture,  in  this  case,  is  associated  with  comminution  (i.e.  formation 
of  a  large  number  of  sub-millimeter  size  fragments).  Nevertheless, 
under  both  low-  and  high-loading  rates,  damage-initiation  as  well 
as  failure/fracture  are  believed  to  be  controlled  by  pre-existing 


material  defects/flaws.  These  flaws  can  become  propagating  cracks 
when  subjected  to  sufficiently  high  stresses.  More  details  regard¬ 
ing  the  basic  physics  of  glass  deformation  and  fracture  behavior 
in  the  quasi-static  and  dynamic  loading  conditions  can  be  found 
in  Refs.  [4,5]. 

2.2.  Key  components  of  the  material  model 

A  detailed  analysis  of  the  key  components  of  the  material  model 
of  Grujicic  et  al.  [4,5]  including  failure  (bulk  and  surface)  probabil¬ 
ity  distribution  function,  crack  nucleation  criteria,  crack  growth 
kinetics,  damage  initiation  and  evolution  equations,  etc.  can  be 
found  in  Refs.  [4,5]. 

2.3.  Mathematical  formulation  of  the  model 


2.3.1.  Coarse-fragmentation  failure/fracture  mode 

As  discussed  earlier,  under  quasi-static  loading  conditions, 
soda-lime  glass  typically  fails  via  the  coarse-fragmentation  mode 
associated  with  the  nucleation  and  propagation  of  a  few  (macro) 
cracks.  This  behavior  of  glass  was  rationalized  as  follows:  Under 
low-loading  rates,  the  rate  of  increase  of  stress  is  also  low.  Hence, 
when  the  first  crack  nucleates  it  begins  to  propagate  (at  a  terminal 
velocity,  typically  equal  to  0.2-0.4  of  the  corresponding  sound 
speed)  and  grow  its  shielding  zone  (within  which  the  stresses  are 
reduced  and  the  potential  crack  nucleating  flaws  made  inactive). 
Meanwhile,  stresses  in  the  unshielded  regions  of  the  material  do 
not  frequently  reach  high-enough  values  to  cause  nucleation  of 
new  cracks.  Consequently,  the  material  fails  via  the  coarse-frag¬ 
mentation  mode. 

Clearly,  in  this  loading  regime  fracture  strength  corresponds  to 
the  stress-level  at  which  the  first  crack  forms.  The  stress  at  which 
the  first  crack  forms,  on  the  other  hand,  is  a  function  of  the  po¬ 
tency/size  of  the  crack-nucleating  flaw.  Consequently,  the  material 
fracture  strength  becomes  a  stochastic  quantity  which  is  defined 
by  an  appropriate  failure-distribution  function,  rather  than  by  a 
mean  value.  Within  the  soda-lime  glass  material  model  of  Grujicic 
et  al.  [4,5],  the  coarse-fragmentation  mode  failure-probability 
function  is  represented  using  a  two-parameter  Weibull-distribu- 
tion  function.  This  function  was  derived  by  treating  the  crack 
nucleation  phenomenon  as  a  Poisson  point-process  and  is  given 
by  the  following  expression: 

PF  =  1  -  exp[AtZ]  (1) 


where  the  failure  probability  PF  represents  the  probability  of  finding 
at  least  one  crack-nucleating  defect  in  the  material  region/domain 
of  size  Z  and  is  the  stress-dependent  defect-density  and  is 
expressed  as: 


(2) 


where  10  is  the  reference  defect  density,  S0  is  a  stress  normalizing 
parameter  and  m  is  the  Weibull  modulus. 

The  failure  probability  distribution  function  (which  is  equal  to 
the  fracture-strength  cumulative  distribution  function)  is  obtained 
by  substituting  Eq.  (2)  into  Eq.  (1 ).  The  fracture  strength  probability 
density  function  is  then  obtained  by  differentiating  the  resulting 
equation  with  respect  to  sigma.  It  should  be  noted  that  despite  an 
explicit  appearance  of  three  parameters  (20,  S0,  m),  the  Weibull  dis¬ 
tribution  function  used  contains  only  two  independent  parameters 
(m  and  S0/2j/m).  Next,  using  the  standard  mathematical  relations 
for  the  single-variant  mean  value  and  standard  deviation,  the 
following  expressions  for  the  fracture-strength  mean  value  and 
standard  deviation  are  derived  [18]: 
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where  T  denotes  the  gamma  function  which  is  defined  as: 

r(x)  =  (x  — 1)!  (5) 

As  discussed  earlier,  to  model  the  growth  of  macro-cracks 
beyond  the  length  of  a  single  element,  the  stress-based  fracture 
initiation  criterion  defined  above  is  combined  with  a  linear-elastic 
fracture-mechanics  based  crack  growth  criterion.  The  procedure 
used  to  implement  the  later  criterion  involves  the  following  steps: 

(a)  first,  adjacent  failed  elements  aligned  in  a  particular  direction 
are  used  to  define  the  associated  crack  length  (2a)  in  that  direction. 
Seven  possible  directions  are  considered,  three  of  which  are 
aligned  with  the  edges  and  the  remaining  four  with  the  diagonals 
of  the  cube-shaped  finite  element;  (b)  for  an  element  located  at  a 
crack  tip,  the  stress  intensity  factor  I<,  is  calculated  by  multiplying 
its  maximum  principal  stress  with  a  factor  y/na;  (c)  when  the  con¬ 
dition  K,  >  Kic,  where  I<lc  is  the  mode  I  critical  stress  intensity  factor 
is  satisfied,  the  fracture  of  the  element  is  initiated.  As  in  the  case  of 
the  stress-based  fracture  criterion,  the  complete  fracture  of  the  ele¬ 
ment  and  its  removal  is  governed  by  the  crack  terminal  velocity. 


2.3.2.  Fine-fragmentation/comminution  failure/fracture  mode 

As  discussed  earlier,  under  dynamic  loading  conditions,  soda- 
lime  glass  typically  fails  via  the  fine-fragmentation  mode  associ¬ 
ated  with  the  nucleation  and  propagation  of  numerous  (micro) 
cracks.  This  behavior  of  glass  was  rationalized  as  follows:  Under 
high-loading  rates,  the  rate  of  increase  of  stress  is  also  high.  Hence, 
before  the  first  nucleated  cracks  can  advance  (and  enlarge  their 
shielding  zones)  considerably,  stresses  in  unshielded  material 
regions  becomes  high  enough  to  cause  nucleation  of  numerous 
cracks.  Consequently,  the  material  fails  via  the  fine-fragmentation 
mode.  It  should  be  also  noted  that  the  neighboring  cracks  with 
compatible  opening  modes  and  orientations  may,  via  their  shield¬ 
ing  zone,  mutually  terminate  each  other’s  growth,  giving  rise  to 
relatively  short  cracks. 

Based  on  the  description  provided  above,  it  was  concluded  that 
failure  and  ultimate  fracture  of  glass  in  this  loading  regime  should 
be  modeled  as  a  progressive-damage  process  rather  than  a  brittle- 
fracture  process  (adopted  in  the  coarse  fragmentation  mode).  This 
was  done  in  the  soda-lime  glass  material  model  of  Grujicic  et  al. 
[4,5].  The  specifics  of  this  material  model  in  the  high-loading  rate 
regime  are  discussed  in  the  remainder  of  this  section  and  involve 
the  following  points: 

(a)  When  a  crack-nucleating  flaw  is  activated  and  the  crack 
begins  to  grow,  a  crack-surrounding  shielding  zone  is 
formed  and  increases  in  size.  Within  the  model,  stresses 
within  the  shielding  zone  are  assumed  to  relax  to  zero 
and,  hence,  activation  of  any  flaw  residing  within  this  zone 
will  be  suppressed. 

(b)  The  size  of  the  shielding-zone  (assumed  to  be  of  a  spherical 
shape  for  bulk  flaws  and  of  a  circular  shape  for  surface 
cracks)  is  postulated  to  increase  during  crack  propagation, 
in  a  self-similar  manner.  Therefore,  the  shielding-zone  size 
at  a  time  f  associated  with  a  crack  which  was  nucleated  at 
the  time  t  is  defined  as: 

Zsh(t,x)  =  S(kC(t-T)]n  (6) 

where  C  =  [E/p]0'5  is  the  speed  of  sound,  £  the  Young’s  modulus,  p 
the  mass  density,  k  =  0.2-0.4,  the  crack  terminal  speed  to  sound 
speed  ratio,  n  is  the  domain-dimensionality  factor  for  the 


fracture-controlling  flaws  (=2,  for  surface  failure  and  =  3,  for  bulk 
failure)  and  S  is  a  shielding  zone  shape  factor  (=47t/3,  for  a  spherical 
bulk  zone  and  n  for  a  circular  shielding  zone); 

(c)  Due  to  the  stress-relaxation  within  the  shielding  zones,  dis¬ 
tinction  is  made  between  the  potentially  crack-nucleating 
(unshielded)  and  the  deactivated  (shielded)  flaws  so  that 
the  total  flaw  density  can  be  decomposed  as: 

A  —  ^non-sh  ^sh  (7) 

where  both  Anon.sh  and  2sj,  are  defined  by  dividing  the  corresponding 
number  of  flaws  by  the  total  domain  size; 


(d)  The  rate  of  loading  affects  the  relative  magnitude  of  the 
unshielded  and  shielded  flaw  densities  as  follows:  (i)  under 
low-loading  rates,  a  large  fraction  of  the  flaws  will  become 
inactive  (due  to  the  presence  of  large  shielding  zone(s)) 
and  under  quasi-static  loading  conditions  all  flaws  (except 
for  the  one(s)  associated  with  the  nucleation  of  first  crack(s)) 
will  become  shielded.  This,  in  turn,  would  yield  a  relatively 
large  value  of  2S/,;  and  (ii)  under  high-loading  rate  condi¬ 
tions,  the  shielding  zone(s)  are  relatively  small  leading  to  a 
relatively  large  value  of  An0n-sh- 

(e)  In  contrast  to  the  coarse-fragmentation  mode  for  which  the 
mean  fracture  strength,  defined  by  Eq.  (3),  is  assumed  to  be 
constant  (but  to  take  different  values  in  the  surface  and  bulk 
material  regions),  the  dynamic  fracture  strength  is  found  to 
be  a  loading-rate  dependent  quantity.  The  functional  form 
for  the  dynamic  fracture  strength  is  given  below. 

(f)  To  derive  an  expression  for  the  dynamic  fracture  strength,  a 
uniform  constant-loading  rate  (&  =  constant)  case  was  ana¬ 
lyzed  first.  Also  distinction  is  made  between  the  externally- 
applied  macroscopic  stress,  Z,  and  its  internal  counterpart 
a  =  at,  where  t  is  the  duration  of  loading.  In  accordance 
with  the  aforementioned  assumption  regarding  zero-stress 
within  the  shielding  zones,  only  non-shielded  portions  of 
the  brittle-material  structure  are  associated  with  the  inter¬ 
nal  stress  level  a.  To  derive  a  relationship  between  the  exter¬ 
nal  stress  (pertains  to  the  entire  bulk  and  surface  domains) 
and  the  internal  stress  (pertains  only  to  the  unshielded  por¬ 
tions  of  these  domains)  a  scalar  damage  parameter,  D,  was 
defined  as  a  ratio  of  the  union  of  all  shielding-zone  vol¬ 
umes/surfaces  and  the  structure  volume/surface.  This 
yielded  the  following  relation: 


Z  =  <7(t)(l  -  D(o))  =  <7t(l  -  D(&,  t)) 


(8) 


where  D  is  taken  to  depend  on  a  and  f,  since  the  product  of  the  two 
is  equal  to  a  and  affects  the  density  of  crack  nucleating  flaws  via  Eq. 
(2),  while  t  alone  affects  the  size  of  the  shielding  zones  via  Eq.  (6). 

According  to  Eq.  (8),  as  the  loading  duration  increases,  the  a(t ) 
increases  causing  Z  to  increase  while  the  (1  —  D(a))  term 
decreases  causing  Z  to  decrease.  Thus,  there  is  a  critical  loading 
duration  at  which  Z  reaches  a  maximum  value  and  this  value, 
2- max*  is  defined  as  the  material  dynamic  strength.  The  material  dy¬ 
namic  strength  is  then  defined  as: 


(9) 


Before  Eq.  (9)  could  be  utilized,  an  expression  is  required  for  the 
internal-stress  dependence  of  the  damage  parameter.  Towards  that 
end  and  following  Denoual  and  Hild  [20,22],  the  damage  parame¬ 
ter  D  is  set  equal  to  the  probability  of  flaw  shielding  (i.e.  the  prob¬ 
ability  of  finding  a  flaw  within  a  shielding  zone),  Psh,  which  is  given 
in  accordance  with  Eq.  (1)  as: 


D=Psh  =  1  -  exp (-AtZsh) 


(10) 
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where  Zsu  is  the  average  size  of  the  shielding  zone  defined  as: 


Mt)zsh(t) 


dlt 

dt 


[kC(t  -  T)]"dT 


(11) 


When  deriving  Eq.  (11),  it  was  taken  into  account  that  the 
shielding  zone  average  size  Zsh  at  time  t,  when  the  total  flaw  den¬ 
sity  is  2t(t),  the  crack  in  question  could  have  been  nucleated  at  any 
time,  t  with  0  ^  t  ^  t  and  that  the  corresponding  shielding-zone 
size  is  [kC(t  -  t)]n.  In  general,  the  probability  for  crack  nucleation 
is  not  constant  in  the  [0,t]  time  interval  but  rather  scales  with 
the  rate  of  activation  of  the  flaws  at  time  t,  given  as  with 

Jo  •Mt)  ^|Tdr  =  1  (since  for  a  shielding  zone  to  exist  the  crack  must 
have  nucleated  at  some  time  in  the  [0,  t]  time  interval). 

In  the  case  of  uniform  loading  under  constant  stress  rate  condi¬ 
tions  considered  here,  and  via  Eq.  (2),  the  term  can  be  written 
as: 


dlt 

dt 


20m&mtm  1 


(12) 


After  substitution  of  Eq.  (12)  into  Eq.  (11)  and,  in  turn,  into  Eq. 
(10),  and  subsequent  integration,  the  following  internal-stress 
dependent  damage-parameter  function  is  obtained: 


D  =  1  -  exp 


( 


\ 


mini 


m+n\ 

Oc  ' 

(m  +  n)! 


where  the  characteristic  internal  stress  ac  is  defined  as: 


Gc  =  Gtc  = 


Sq  dn 
UoS(kC)" 


(13) 


(14) 


Substitution  of  Eq.  (13)  into  Eq.  (8)  and  differentiation  of  the 
resulting  equation  in  accordance  with  Eq.  (9)  yields: 


<7 


=  (JS 


m  (Jn 


\X0S(kC)n 

and 

G [  dynamic  ~  Zmax  —  Gc 


(m  +  n-  1)! 
mini 


1  (m  +  n-  1)! 
e 


mini 


(15) 


(16) 


where  oy dynamic  is  the  material  average  fracture  strength  in  the  fine- 
fragmentation  mode.  Substitution  of  Eq.  (15)  into  Eq.(13)  yields  the 
following  expression  for  the  extent  of  material  damage  at  the  onset 
of  cracking  under  high-loading  rates: 


D|v  =1- 

lz-max 


\  m+n 

ej 


(17) 


Furthermore,  setting  a  =  at  in  Eq.  (15)  gives  the  duration  of  loading 
at  the  onset  of  dynamic  fracture  as: 


/  \m+"  /(m  +  n  -  l)!\”+n 

t|rm“  =  V;.0S(kC)"<7m )  V  S  ) 


(18) 


2.3.3.  Comparison  between  quasi-static  and  dynamic  failure  initiation 
To  help  understand  the  relationship  between  coarse-fragmenta- 
tion  and  fine-fragmentation  failure/fracture  modes  for  soda-lime 
glass,  fracture  strength  vs.  constant  stress-rate  plot  is  constructed 
and  displayed  in  Fig.  la  and  b  for  the  bulk  and  surface  failure, 
respectively.  To  assist  in  interpretation  of  these  figures,  a  second 
horizontal  axis,  Zejj/Zc,  is  introduced  where  Zc  is  the  Z  equivalent 
of  ac.  In  both  Fig.  la  and  b,  two  fracture  strength  curves  are 
displayed:  (a)  one  corresponding  to  the  average  stress-rate  inde¬ 
pendent  quasi-static  fracture  strength,  Eq.  (3);  and  (b)  the  average 
dynamic  fracture  strength  (increases  with  an  increase  in  the  stress- 


ZJZC,  No  Units 


Z  /Z  ,  No  Units 

eff  C 


Fig.  1.  The  transition  between  the  quasi-static  (coarse  fragmentation)  and  the 
dynamic  (fine  fragmentation)  brittle-fracture  modes  as  a  function  of  an  increase  in 
stress  rate  for:  (a)  bulk;  and  (b)  surface  mode  of  fracture.  See  text  or  the  explanation 
of  Zeff  and  Zc. 


rate),  Eq.  (16).  The  relevant  quasi-static  and  dynamic  facture  mate- 
rial-model  parameters  used  in  the  construction  of  Fig.  la  and  b  are 
listed  in  Table  1. 

It  should  be  noted  that  in  Fig.  la  and  b,  the  two  expressions  for 
the  (mean)  fracture  strength,  Eqs.  (3)  and  (16),  are  valid  only  over  a 
limited  range  of  stress  rates  and  that  the  ranges  are  different  for 
the  two  relations.  That  is,  in  the  high  stress  rate  range,  defect 
shielding  is  less  pronounced  and,  hence,  the  coarse-fragmentation 
fracture  strength  relation,  Eq.  (3),  (which  was  based  on  an  assump¬ 
tion  of  complete  shielding  of  the  surrounding  material  by  the  first 
nucleated  crack)  is  not  valid.  On  the  other  hand,  the  dynamic  frac¬ 
ture  strength  relation,  Eq.  (16),  is  not  valid  in  the  low  stress  rate 
range  (i.e.  at  lower  values  of  Zejj/Zc),  since  in  this  case  the  size  of 
the  accompanying  shielding  zone  at  the  onset  of  fracture  exceeds 
the  total  domain  size.  The  approximate  stress-rate  ranges  of  valid- 
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Table  1 

Mechanical  property  parameters  for  soda-lime  glass  used  in  the  present  work. 


Property 

Symbol 

Value 

Unit 

Young’s  modulus 

E 

70.0 

GPa 

Poisson’s  ratio 

V 

0.22 

N/A 

Density 

p 

2500 

kg/m3 

Mean  fracture  toughness 

K,c 

0.75 

MPa  m1'2 

Weibull  parameters  for  surface  controlled  fracture 

Weibull  modulus 

m 

7 

N/A 

Mean  static  fracture  strength 

G /static 

100 

MPa 

Effective  surface 

0.01 

m2 

Weibull  scale  parameter3 

So/4/m 

55.3 

MPa  m2/m 

Weibull  parameters  for  volume  controlled  fracture 

Weibull  modulus 

m 

30 

N/A 

Mean  static  fracture  strength 

G /static 

230 

MPa 

Effective  volume 

to-4 

m3 

Weibull  scale  parameter3 

c  /  il/m 

So/V 

186 

MPa  m2,m 

a  Computed  using  Eq.  (3). 


ity  for  the  two  fracture-strength  functions  are  indicated  using 
“solid"  lines  in  Fig.  la  and  b.  The  intermediate  stress-rate  range 
is  represented  by  the  “dashed"  portions  of  the  curves.  In  general, 
within  the  intermediate  stress-rate  range,  one  would  expect  a 
monotonic  transition  between  the  low  and  high  stress-rate 
branches.  However,  in  the  model  of  Grujicic  et  al.  [4,5]  the  transi¬ 
tion  between  the  coarse-fragmentation  mode  and  the  fine-frag¬ 
mentation  mode  was  assumed  to  occur  abruptly  at  a  constant 
value  of  the  stress  rate  (defined  as  the  point  of  intersection  of 
the  two  fracture-strength  curves  in  Fig.  la  and  b). 


2.3.4.  Damage-induced  material  property  degradation 

As  discussed  earlier,  fine-fragmentation  mode-induced  damage 
causes  degradation  of  the  material  stiffness  and  strength.  In  order 
to  quantify  the  progression  of  material-property  degradation,  a 
damage  evolution  equation  is  required.  Such  an  equation  is 
obtained  by  differentiating  Eq.  (13)  with  respect  to  a,  to  get: 


dD/da 


m\n\(m  +  n)crm+n~1 

(m  +  n)!er™+n  (  "  > 


(19) 


It  should  be  noted  that,  within  the  soda-lime  glass  material 
model  of  Grujicic  et  al.  [4,5],  damage  is  assumed  to  be  of  an  isotro¬ 
pic  character  so  that  the  material  remains  isotropic  during  failure. 
Out  of  the  two  isotropic  material  elastic  constants  (Young’s  modu¬ 
lus  and  Poisson’s  ratio),  only  the  Young’s  modulus  is  assumed  to 
degrade. 

Young’s  modulus  of  glass  is  then  degraded  according  to  the 
following  relation: 


E  =  E0(l  -  D) 


(20) 


where  subscript  0  is  used  to  denote  a  quantity  pertaining  to  glass  in 
its  virgin  state. 

A  comparison  of  Eqs.  (8)  and  (20)  shows  that,  in  the  multiple- 
fragmentation  mode,  the  unshielded  portions  of  the  material  are 
associated  with  the  internal  stress  er  and  the  virgin-material  Young’s 
modulus,  E0,  while  the  entire  unshielded  +  shielded  material  region 
is  associated  with  the  macroscopic  stress,  I  =  <r(l  -  D(cr))  and  the 
degraded  Young's  modulus,  E  =  E0(t  -  D). 


2.4.  Material-model  parameterization 

A  summary  of  the  soda-lime  glass  material-model  parameters 
used  in  the  present  work  is  provided  in  Table  1.  This  summary  con¬ 
tains  the  parameters  assessed  by  Grujicic  et  al.  [4,5]  using  the 
open-literature  experimental  and  computational  results  and  newly 


Fig.  2.  Weibull-type  failure  probability  curves  for  soda-lime  glass  in  the  case  of 
bulk  and  surface  modes  of  fracture. 


assessed  parameters.  The  new  assessment  involved  the  same 
open-literature  data,  but  it  was  of  a  higher  rigor. 

The  virgin-material  elastic  stiffness  properties  (including 
Young's  modulus  and  Poisson’s  ratio),  density,  and  the  mode-I  crit¬ 
ical  stress  intensity  factor  listed  are  taken  from  Grujicic  et  al.  [4,5]. 
The  same  is  true  for  the  bulk-controlled  fracture  parameters.  The 
surface-controlled  fracture  parameters  are  obtained  in  the  new 
assessment  procedure.  The  failure  probability  functions  for  the 
bulk  and  the  surface  modes  of  failure  based  on  the  parameteriza¬ 
tion  provided  in  Table  1  and  a  f  cm3  (cube)  soda-lime  specimen 
are  displayed  in  Fig.  2.  As  expected,  surface  fracture  is  associated 
with  a  lower  level  of  the  fracture  strength  and  a  broader  distribu¬ 
tion.  As  stated  earlier,  transition  between  the  coarse-fragmentation 
and  fine-fragmentation  brittle-fracture  modes  is  assumed  to  take 
place  at  a  constant  stress-rate.  Based  on  the  results  displayed  in 
Fig.  la  and  b,  this  stress-rate  was  set  to  a  value  of  2.68  MPa/ps 
for  the  bulk  facture  and  1.83  MPa/ps  for  the  surface  fracture. 

In  addition  to  the  parameters  listed  in  Table  1,  the  following 
conditions  were  defined:  (a)  in  accordance  with  the  previous  dis¬ 
cussion,  the  crack  terminal  velocity  is  set  to  a  value  equal  to  30% 
of  the  material  sound  speed  (the  latter  is  calculated  using  the 
materials  Young's  modulus  and  density);  and  (b)  as  stated  earlier, 
depending  on  the  (bulk  vs.  surface)  location  of  the  flaws,  the  crack 
shielding  zones  are  assumed  to  be  either  spherical  or  circular, 
respectively.  Consequently,  for  the  two  cases,  the  dimensionality 
factor,  n,  and  the  shielding  zone  shape  factor,  S,  are  set  to  3  and 
4/37t  for  bulk  failure,  and  2  and  7t  for  surface  failure. 

As  far  as  the  linear-elastic  fracture-mechanics  based  macro¬ 
crack  growth  model  is  concerned,  it  is  associated  with  a  single 
material  parameter,  i.e.  the  mode-I  critical  stress  intensity  factor, 
KIC.  In  accordance  with  the  stress-based  macro-cracking  initiation 
criterion,  the  mode-I  critical  stress  intensity  factor  was  taken,  in 
the  prior  work  of  Grujicic  et  al.  [4,5],  to  be  a  stochastic  quantity 
and  given  by  a  Weibull-distribution  function  with  the  same  refer¬ 
ence  defect  density,  X0,  and  Weibull  modulus,  m,  (have  different 
values  in  the  bulk  and  surface  regions)  as  in  the  fracture  strength 
case.  As  far  as  the  scaling  parameter  is  concerned,  it  was  computed 
using  an  equation  analogous  to  Eq.  (3)  and  assuming  a  constant 
mean  value  of  the  mode-I  critical  stress  intensity  factor  of 
0.75  MPa  m1/2  [20]  (for  both  the  bulk  and  surface  regions).  In  the 
present  work,  I<IC  was  not  considered  to  be  a  stochastic  quantity. 
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This  is  not  considered  to  be  a  simplification,  but  rather  a  physically 
more  correct  treatment  of  this  quantity.  That  is,  in  brittle  solids 
such  as  glass,  1<IC  is  related  to  the  material  surface  energy  and  the 
latter  quantity  (in  the  absence  of  large  microstructure/composi¬ 
tional  heterogeneities)  is  generally  considered  to  be  deterministic. 

3.  Soda-lime  glass  non-linearity  and  inelasticity 

As  explained  earlier,  the  main  objective  of  the  present  work  is  to 
examine  (computationally)  the  role  of  material  non-linearity  and 
inelastic  deformation  behavior  on  the  dynamic  response  of 
soda-lime  glass.  In  this  section,  a  brief  description  is  provided  of 
the  basic  physics  for  the  two  phenomena,  their  mathematical 
formulation,  and  their  implementation  into  the  soda-lime  material 
model  of  Grujicic  et  al.  [4,5]. 

3.1.  Underlying  physics 

In  the  original  model  of  Grujicic  et  al.  [4,5],  the  material  volu¬ 
metric  response  was  assumed  to  be  linear  elastic,  both  in  compres¬ 
sion  and  tension.  This  assumption  is  not  fully  consistent  with  a 
number  of  experimental  findings  [14,15]  which  showed  that  under 
sufficiently  high  pressures  (greater  than  ca.  3.5  GPa)  the  compres¬ 
sive  loading  response  of  soda-lime  glass  becomes  visibly  nonlinear 
(material  nonlinearity)  and  irreversible  (inelastic).  This  behavior  of 
soda-lime  glass  is  depicted  in  the  axial  stress  vs.  normalized  spe¬ 
cific  volume  plot,  Fig.  3,  obtained  in  a  standard  flyer  plate  experi¬ 
ment  [28].  Under  hydrostatic  tension,  however,  material  fails  at  a 
relatively  low  pressure  (less  than  ca.  300  MPa)  and,  prior  to  the 
onset  of  failure,  the  material  response  is  linear  elastic. 

A  literature  review  carried  out  as  part  of  the  present  work 
established  that  there  is  no  general  consensus  regarding  the  nature 
of  the  material  intrinsic  phenomena  and  processes  responsible  for 
the  aforementioned  nonlinearity  and  inelasticity  of  soda-lime 
glass.  The  following  main  theories  appear  noteworthy:  (a)  shear  in¬ 
duced  micro-cracking  [29].  The  main  supporting  evidence  for  this 
phenomenon  is  the  experimental  observation  of  reduced  spallation 
strength  in  (compressively)  pre-shocked  soda-lime  glass;  (b)  irre¬ 
versible  compaction  characterized  by  permanent  density  increase 
and  non-significant  changes  in  the  material  bond  structure  and 
molecular-level  topology  e.g.  [30].  The  main  supporting  evidence 
in  this  case  comes  from  post-mortem  characterization  of  the  mate- 


Fig.  3.  Typical  axial-stress  vs.  normalized  specific  volume  response  during 
dynamic-compressive  loading  and  unloading. 


rial  density  and  microstructure;  and  (c)  various  phase  transforma¬ 
tions  which  are  accompanied  with  both  significant  irreversible 
increases  in  material  density  and  major  changes  in  material  bond 
structure  and  molecular-level  topology  [16].  As  in  the  previous 
case,  the  main  supporting  evidence  comes  from  post-mortem  char¬ 
acterization  of  the  material  density  and  microstructure.  In  addi¬ 
tion,  molecular-level  simulations  of  shock  wave  generation  and 
propagation  in  soda-lime  glass  provided  evidence  for  changes  in 
the  Si-atom  coordination  and  Si-0  ring  topology  in  the  as-shocked 
material  state.  An  example  of  such  molecular-level  computational 
results  obtained  in  our  ongoing  work  is  displayed  in  Fig.  4a  and  b. 

3.2.  Mathematical  formulation 

As  discussed  above,  the  origin  of  soda-lime  glass  non-linearity 
and  inelasticity  has  not  yet  been  fully  established.  This,  however, 


Fig.  4.  Changes  in  the  soda-lime  glass  bond  structure  and  molecular-level  topology 
under  shock  loading:  (a)  development  of  fivefold  coordinated  silicon  atoms 
(highlighted  in  green);  and  (b)  formation  of  smaller  Si-0  rings  (highlighted  in 
purple  and  yellow).  (For  interpretation  of  the  references  to  color  in  this  figure 
legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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is  not  a  hindrance  from  the  point  of  implementation  of  these  phe¬ 
nomena  as  an  enhancement  to  the  continuum-level  material  mod¬ 
el  for  soda-lime  glass  developed  by  Grujicic  et  al.  [4,5].  In  other 
words,  by  examining  the  experimental  data  published  in  the  open 
literature,  it  was  possible,  within  the  present  work,  to  quantify  the 
extent  of  material  nonlinearity  and  deformation  irreversibility  as  a 
function  of  the  maximum  compressive  pressure.  It  should  be  noted 
that  the  shear  component  of  the  material  response  is  assumed  to 
remain  linear.  In  the  remainder  of  this  section,  brief  explanations 
are  provided  of  the  mathematical  formulations  used  to  quantify 
soda-lime  glass  nonlinearity  and  inelasticity. 

The  experimental  data  published  in  the  open  literature  [28]  are 
used  to  determine  the  dependence  of  soda-lime  glass  pressure  on 
the  normalized  specific  volume  in  the  form: 

(  k,  +k2{v/v0)  +  k3(v/v 0)2,  for  v/v0  >  0.916 
P( v/v0)  =  l  k4,  for  0.813  <  v/v0  <  0.916  (21) 

l  /<5  +  k6(v/v 0)  +  k7(v/v0)2,  for  v/v0  <  0.813 

where  k,  =  -301.3  GPa,  k2  =  668.9  GPa,  k3  =  - 367.6  GPa,  k4  = 
2.97  GPa,  k5  =  780  GPa,  k6  =  -1920  GPa  and  k7  =  1184  GPa.  Exami¬ 
nation  of  Eq.  (21)  reveals  the  presence  of  three  distinct  v/v0 
(compressive)  regions:  (a)  0.916  <  vlv0<  1.0,  the  initial  anomalous 
compression  regime  characterized  by  a  continuous  volumetric  soft¬ 
ening;  (b)  0.813  <  v!v0<  0.916,  the  irreversible  compaction  regime; 
and  (c)  v/v0  <  0.813,  the  normal  compression  regime  characterized 
by  continuous  volumetric  stiffening.  Coefficients  kj  to  k7  are 
obtained  by  applying  the  conventional  least-squares  curve-fitting 
procedure  to  the  experimental  data  reported  in  Ref.  [28]. 

Eq.  (21 )  is  used  in  the  present  work  in  two  different  ways:  (a)  to 
define  volumetric  loading  and  unloading  paths  under  the  assump¬ 
tion  that  the  material  is  non-linear.  In  this  case,  the  tangent  bulk 
modulus  is  defined  using  the  negative  slope  of  the  P  vs.  v/v0  curve; 
and  (b)  to  define  the  loading  path  under  the  assumption  that  the 
material  is  linear-elastic  and  volumetrically-plastic.  In  other 
words,  the  P  vs.  v/v0  curve  represents  the  irreversible  volumetric- 
strain  hardening  behavior  of  soda-lime  glass.  Unloading  and  elas¬ 
tic  reloading,  on  the  other  hand,  are  taken  to  be  associated  with  a 
linear-elastic  material  response  and  a  volumetric-strain  indepen¬ 
dent  bulk  modulus  (equal  to  the  bulk  modulus  of  the  virgin  mate¬ 
rial,  Table  1). 

3.3.  Implementation  into  a  user-material  subroutine 

Mathematical  formulation  for  the  material  nonlinearity  and 
deformation  inelasticity  are  used  to  enhance  the  continuum-level 
material  model  for  soda-lime  glass  [4,5],  The  resulting  enhanced 
material  model  was  coded  using  the  Intel  Fortran  computational 
language  and  implemented  as  a  VUMAT  Material  User  Subroutine 
within  the  commercial  finite  element  program  ABAQUS/Explicit 
[31].  At  the  beginning  of  each  analysis,  the  material  subroutine  is 
compiled  and  linked  with  the  finite  element  solver  which  enables 
ABAQUS/Explicit  to  pass  the  necessary  material  state  and  deforma¬ 
tion  variables  to  the  material  model  for  each  element  integration 
point  at  each  time  step  and  subsequently  received  back  the  up¬ 
dated  material  state  and  stress  variables. 

The  basic  procedure  for  coupling  the  ABAQUS/Explicit  finite- 
element  solver  with  the  VUMAT  Material  User  Subroutine  at  each 
time  increment  at  each  integration  point  of  each  element  can  be 
outlined  as  follows: 

(a)  The  finite-element  integration-point-specific  previous  time- 
increment  stresses  and  material  state  variables  as  well  as  the 
current  time-step  incremental  strains  are  passed  to  the 
VUMAT  by  the  ABAQUS/Explicit  finite-element  solver.  Spe¬ 
cifically,  the  glass  material  model  formulation  used  in  the 


present  work  requires  the  following  state  variables  to  be 
passed  into  the  VUMAT  in  addition  to  the  requisite  strain 
increment  components:  (i)  material  initial  strength  value, 
assigned  randomly  from  the  material-specific  Weibull  distri¬ 
bution  (although  this  quantity  does  not  change,  it  is  treated 
as  a  state  variable  for  practical/operational  reasons);  (ii)  the 
extent  of  coherent  damage;  (iii)  minimum  specific  volume 
achieved  (to  track  irreversible  densification);  and  (iv)  a 
variable  defining  the  deletion  status  of  the  element;  and 
(b)  Using  the  material  state  information  passed  to  the  VUMAT  in 
(a),  and  the  soda-lime  glass  material  model  presented  in 
Section  2,  the  material  stress  state  as  well  as  the  updated 
values  of  the  material  state  variables  are  calculated  and 
returned  to  the  ABAQUS/Explicit  finite-element  solver.  Addi¬ 
tional  information  that  is  required  to  be  updated,  but  not 
passed  back  to  the  solver,  in  order  to  handle  the  macro¬ 
cracking  component  of  the  material  modeling  includes  the 
element  location,  failure  status,  and  the  direction  of  crack 
opening  (where  appropriate)  which  are  stored  in  the  form 
of  globally  available  matrices. 

4.  Problem  formulation  and  computational  analysis 

4.1.  Problem  formulation 

To  examine  the  effect  of  material  non-linearity  and  inelasticity 
on  the  dynamic  behavior  of  soda-lime  glass,  the  case  of  a  symmet¬ 
ric  flyer  plate  impact  is  analyzed.  Within  this  problem,  a  plate-like 
target  is  impacted  by  a  plate-like  projectile  at  a  zero  obliquity  an¬ 
gle  (normal  impact).  The  target  and  projectile  are  composed  of  the 
same  (soda-lime  glass)  material.  The  projectile  thickness  (5  mm,  in 
the  present  case)  is  selected  to  be  less  than  the  target  thickness 
(10  mm,  in  the  present  case)  to  ensure  that  the  so-called  candidate 
spallation  plane  lies  within  the  target.  In  the  limit  of  a  large  ratio  of 
the  target/projectile  lateral  dimensions  to  the  thickness,  the  prob¬ 
lem  can  be  treated  as  being  one-dimensional  in  character. 

A  schematic  of  the  typical  time,  t  vs.  (Lagrangian-type)  spatial 
coordinate,  X,  plot  for  the  problem  at  hand  is  depicted  in  Fig.  5. 
Examination  of  this  figure  reveals  the  following:  (a)  at  the  moment 
of  impact  two  centered  simple  waves  are  generated  at  the  impact 
surface  and  propagate  toward  the  projectile  and  target  free  sur¬ 
faces.  The  reason  that  the  impact  does  not  generate  shock  waves 
is  that  soda-lime  glass  shows  an  anomalous  concave  downward 
P  vs.  v/v0  behavior  (at  least  within  a  0-3  GPa  range,  analyzed  in 
the  present  work),  Fig.  3;  (b)  upon  the  reflection  of  these  waves 
from  the  target/projectile  free  surfaces,  two  simple  waves  (with 
converging  characteristics)  are  formed;  and  (c)  two  approaching 
simple  waves  intersect  within  a  region  marked  by  vertices  a,  b,  c, 
and  d,  Fig.  5.  The  b,  c,  and  d  bounded  portion  of  this  region  contains 
the  material  which  first  experiences  tensile  stresses.  The  candidate 
spallation  plane  is  defined  as  the  X-location  within  this  region  at 
which  the  tensile  stress  first  becomes  equal  to  the  material  fracture 
strength.  Typically,  spallation  fracture  initiates  on  this  plane. 

4.2.  Computational  analysis 

4.2.1.  Computational/geometrical  domain(s) 

As  explained  above,  the  problem  at  hand  is  of  a  one-dimensional 
character.  The  computational  domain  employed  is  composed  of  two 
Lagrangian  regions,  one  representing  the  projectile  and  the  other 
representing  the  target,  Fig.  6.  The  dimension  of  the  two  regions  in 
the  impact  direction  is  set  equal  to  the  respective  plate  thicknesses 
(defined  above).  The  lateral  dimensions  are  inconsequential  to  the 
computational  results  and  thus  are  assigned  arbitrarily  large  values 
to  facilitate  post-processing  visualization. 
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Fig.  5.  Schematic  of  a  typical  time  vs.  (Lagrangian)  spatial  coordinate  plot  for  a 
symmetric  soda-lime  glass  flyer-plate  impact.  Note,  T  and  L  are  used  to  denote  the 
leading  and  trailing  portions  of  the  wave,  respectively. 

4.2.2.  Finite  element  mesh 

The  two  regions  are  meshed  to  each  form  a  column  of  single  ele¬ 
ments  running  in  the  direction  of  impact.  To  determine  the  axial 
mesh  dimension,  the  mesh  was  gradually  refined  until  further 
refinement  did  not  considerably  affect  the  results.  This  procedure 
yielded  the  meshed  axial  dimension  of  10  pm. 

4.2.3.  Material  model 

Due  to  the  symmetric  nature  of  the  flyer-plate  impact  problem, 
the  projectile  and  the  target  are  composed  of  identical  (soda-lime 
glass)  material.  This  material  is  modeled  represented  using  three 
different  models:  (a)  the  original  model  of  Grujicic  et  al.  [4,5];  (b) 
the  present  rendition  of  the  this  model  which  includes  only  addi¬ 
tional  nonlinear  elastic  effects;  and  (c)  the  present  rendition  of 
the  this  model  which  includes  both  additional  nonlinear  elastic 
and  inelastic  effects. 

4.2.4.  Initial  conditions 

At  zero  simulation  time,  the  projectile  is  assigned  an  incident 
velocity  in  a  100-900  m/s  range,  while  the  target  is  assumed  to 
be  stationary. 

4.2.5.  Boundary  conditions 

Zero  lateral  velocity  boundary  conditions  are  applied  to  all  the 
nodes  in  the  model.  This  boundary  condition  ensures  the  uniaxial 
strain  material  state  that  is  expected  under  flyer-plate  impact 
conditions. 
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Fig.  6.  The  computational  Lagrangian  domain  representing  the  projectile  and  the 
target  plates. 


4.2.6.  Contact  conditions 

The  projectile/target  interaction  is  modeled  using  a  ‘‘penalty’’ 
contact  method  within  which  the  penetration  of  the  surfaces  into 
each  other  is  resisted  by  linear  spring  forces/contact-pressures 
with  values  proportional  to  the  depth  of  penetration.  Within  this 
algorithm,  pressure  is  transmitted  only  when  the  two  bodies  are 
in  contact  and  there  is  no  limit  to  the  magnitude  of  this  contact 
pressure. 

4.2.7.  Computational  method 

The  flyer-plate  impact  problem  described  above  is  analyzed 
computationally  in  the  present  work  using  a  transient,  non-linear 
dynamics,  explicit,  Lagrangian,  finite  element  analysis.  Further 
details  regarding  this  method  can  be  in  our  prior  work  [26],  It 
should  be  noted  that  no  variable  mass  scaling  algorithm  was  used 
to  improve  the  computational  efficiency  due  to  the  relatively  low 
computational  cost.  Furthermore,  due  to  the  very  fine  nature  of 
the  mesh  and  the  absence  of  shock  waves,  no  bulk  viscosity 
algorithm  (aimed  at  mitigating  the  computational  challenges  asso¬ 
ciated  with  large  field  gradients)  was  used. 

5.  Results  and  discussion 

In  this  section,  the  main  results  obtained  in  the  aforementioned 
transient  nonlinear  dynamics  computational  analysis  of  the  sym¬ 
metric  flyer-plate  impact  scenario  are  presented  and  discussed. 
Specifically,  the  results  pertaining  to  the;  (a)  wave  structure/mo¬ 
tion/interaction;  (b)  target  back-face  velocity  history  analysis; 
and  (c)  impact  velocity  spallation  threshold  analysis. 

5.3.  Wave  structure/motion/interaction 

A  typical  particle-velocity  field  plot  (at  8  equally-spaced,  by 
120  ns,  post-impact  times)  for  the  case  of  the  projectile/target 
material  represented  using  the  original  soda-lime  glass  material 
model  of  Grujicic  et  al.  [4,5]  is  displayed  in  Fig.  7.  Within  the 
impact  scenario  depicted  in  Fig.  7  (as  well  as  in  Figs.  8  and  9), 
the  projectile  was  assigned  an  incident  particle  velocity  of 
200  m/s.  Examination  of  Fig.  7  reveals:  (a)  the  formation  of  two 
diverging  compression  waves  emanating  from  the  impact  inter¬ 
face.  The  two  waves  bound  a  projectile/target  region  characterized 
by  a  particle  velocity  of  ca.  100  m/s;  (b)  upon  reflection  of  these 
waves  from  their  respective  free  surfaces,  two  approaching  release 
waves  are  formed.  The  right  propagating  release  wave  leaves  an 
effectively  stationary  material  in  its  wake,  while  the  material  be¬ 
hind  the  left  propagating  release  wave  acquires  the  particle  veloc¬ 
ity  matching  the  projectile  initial  velocity;  (c)  upon  the 
intersection  of  these  two  release  waves,  a  tensile  region  bounded 
by  the  two  (now  diverging)  release  waves  is  formed.  The  existence 
of  tensile  stresses  in  this  region  is  related  to  the  fact  that  the  mate¬ 
rial  particles  at  the  left  boundary  of  this  region  are  stationary  while 
those  at  the  right  boundary  move  (to  the  right)  at  the  projectile  ini¬ 
tial  velocity;  and  (d)  spall  fracture  takes  place  in  this  region  along 
the  candidate  spall  plane(s)  resulting  in  the  formation  of  a  princi¬ 
pal  fragment  (with  a  linear  momentum  nearly  equal  to  that  of  the 
projectile  prior  to  impact)  of  target  material  and  several  smaller 
fragments.  Small  fragments  are  formed  as  a  result  stress  “ringing" 
following  the  formation  of  the  principal  fragment.  While  this  arti¬ 
fact  could  be  prevented  by  employing  an  artificial  viscosity  algo¬ 
rithm,  this  was  not  done  in  the  present  work  in  order  to  preserve 
the  natural  wave  profiles.  While  the  observations  made  are  gener¬ 
ally  consistent  with  t  vs.  X  plot  shown  in  Fig.  5,  there  are  a  few 
points  of  disagreement.  As  will  be  more  clearly  shown  in  Fig.  9, 
these  disagreements  are  mainly  related  to  wave  structure  (wave 
front  profile).  Specifically,  spreading  of  the  compression  waves 
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Fig.  7.  A  time-series  (120  ns  intervals)  of  particle-velocity  field  plots  for  the  case  of 
the  projectile/target  material  represented  using  the  original  soda-lime  glass 
material  model  of  Grujicic  et  al.  [4,5]  and  an  initial  flyer-plate  velocity  of  200  m/s. 


Fig.  8.  A  time-series  (120  ns  intervals)  of  particle-velocity  field  plots  for  the  case  of 
the  projectile/target  material  represented  using  the  nonlinear  elastic  material 
model  and  an  initial  flyer-plate  velocity  of  200  m/s. 

and  steepening  of  the  release  waves  (expected  in  anomalous 
material  like  soda-lime  glass)  is  not  observed  (since  within  the 
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Fig.  9.  Temporal  evolution  of  the  axial-stress  at  three  different  locations  within  the 
target:  (a)  a  point  close  to  the  flyer-plate/target  impact  interface;  (b)  a  point  on  the 
spall  plane;  and  (c)  a  point  close  to  the  target  back-face. 


attendant  linear  elastic  framework,  both  of  these  wave  types  are 
steady,  weak  elastic  waves.). 

The  particle  velocity  field  plot  analogous  to  that  displayed  in 
Fig.  7,  but  for  the  case  of  the  glass  material  model  with  the  afore¬ 
mentioned  nonlinear  elastic  modification,  is  shown  in  Fig.  8.  The 
same  four  general  observations  regarding  the  wave  propagation/ 
reflection/interaction  made  are  mirrored  in  Fig.  8.  On  the  other 
hand,  noticeable  differences  exist  between  the  two  sets  of  results 
regarding  the  wave  structure/profile.  Specifically,  the  results  dis¬ 
played  in  Fig.  8  provide  clear  evidence  of  compression  wave 
spreading  and  release  wave  steepening.  As  mentioned  earlier, 
these  effects  will  be  shown  more  clearly  in  Fig.  9. 

The  corresponding  particle-velocity  field  plots  for  the  case  of 
the  glass  material  model  with  the  aforementioned  nonlinear  elastic 
and  inelastic  modifications  are  not  shown  since  they  are  quite  sim¬ 
ilar  to  the  ones  displayed  in  Fig.  8. 

Temporal  evolution  of  the  axial-stress  at  three  different  locations 
within  the  target  is  displayed  in  Fig.  9a-c.  The  results  displayed  in 
Fig.  9a  pertain  to  an  element  located  near  the  projectile/target  inter¬ 
face,  those  in  Fig.  9b  to  an  element  within  the  spall  region,  while 
those  in  Fig.  9c  pertain  to  an  element  near  the  target  back-face.  In 
each  case,  two  curves  are  displayed  corresponding  to  the  material 
models  employed  in  Figs.  7  and  8.  Comparison  of  the  results  dis¬ 
played  in  Fig.  9a-c  clearly  reveals  spreading  of  the  right  propagating 
compression  wave  for  the  case  of  the  nonlinear  elastic  material 
model.  Furthermore,  examination  of  Fig.  9a  shows  that  the  right 
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propagating  release  wave  has  steepened  considerably  for  the  same 
material  model  case.  No  similar  changes  in  the  wave  profile  are 
apparent  for  the  original  soda-lime  glass  material  model  [4,5].  In 
passing,  it  should  be  noticed  that  upon  the  intersection  of  the  two 
release  waves,  Fig.  9b,  the  stress  first  becomes  tensile,  upon  spalla¬ 
tion,  retreats  to  zero  stress. 

5.2.  Target  back-face  velocity  histoiy 

The  most  common  way  of  determining  experimentally  the 
shock  hugoniot  relations  for  a  material  is  to  employ  flyer-plate 
experiments  and  measure  the  shock  speed  and  the  (upstream) 
particle  velocity  at  different  flyer-plate  impact  velocities.  In  these 
experimental  procedures,  free  surface  particle  velocities  are 
typically  measured  at  the  target  back-face  using  various  laser 
interferometry  techniques.  Subsequently,  a  computational  proce¬ 
dure  is  employed  to  extract  upstream  particle  velocities  from  the 
measured  free-surface  velocities.  In  this  section,  the  results  per¬ 
taining  to  the  history  of  target  back-face  free-surface  velocity  are 
presented  and  discussed. 

Temporal  evolution  of  the  target  back-face  particle  velocities  for 
the  200  m/s  and  900  m/s  projectile  initial  velocity  cases  is 
displayed  in  Fig.  10a  and  b,  respectively.  Three  curves  are  shown 
in  each  figure  corresponding  to  the  three  soda-lime  glass  material 
model  cases. 

Examination  of  the  results  displayed  in  Fig.  10a  and  b  show 
that:  (a)  arrival  of  the  right-propagating  release  wave,  formed  at 
the  spallation-created  free  surface,  to  the  target  (now  the  principal 
fragment)  back-face  produces  the  so-called  "pullback  signal’f  arrest 
and  subsequent  recovery  of  the  back-face  particle  velocity);  (b)  in 
the  lower  flyer-plate  impact  velocity  case,  Fig.  10a,  notable  differ¬ 
ences  are  observed  in  the  free  surface  velocity  history  relative  to 
the  magnitude  of  the  pullback  signal.  However,  no  similar  differ¬ 
ences  are  observed  in  the  ultimate  particle  velocity  achieved  at 
the  target  back-face.  In  the  context  of  the  ballistic  threat  imposed 
by  the  principal  “ flying ”  spall  fragment,  the  magnitude  of  the  pull¬ 
back  is  less  relevant  than  the  ultimate  particle  velocity  attained.  It 
should  be  noted  that  only  two  curves  in  Fig.  10a  are  apparent  due 
to  identical  low  impact  speed  material  response  for  the  nonlinear- 
elastic  and  nonlinear-elastic/inelastic  material  model  cases;  and  (c) 
in  the  higher  flyer-plate  impact  velocity  case,  Fig.  10b,  it  is 
observed  that  all  three  material  model  cases  again  display  a 
pullback  signal,  however,  these  signals  are  quite  comparable  in 
magnitude.  It  should  be  also  noted  that  formation  of  additional 
finer  non-principal  fragments  produces  trailing  pullback  signals, 
in  the  case  of  the  original  soda-lime  glass  material  model  [4,5], 
As  far  as  the  back-face  ultimate  particle  velocities  are  concerned, 
they  are  significantly  higher  in  original  material  model  case. 

5.3.  Impact  velocity  spallation  threshold 

Spallation  is  a  prominent  damage/fracture  phenomenon  which 
is  undesirable  since  it  negatively  affects  the  protection  perfor¬ 
mance  of  a  glass-based  transparent  armor  system  (e.g.  vehicle/ 
structure  interior  fragment  release,  degradation/loss  of  optical 
transparency,  and  decreased  multi-hit  performance).  Hence,  it  is 
important  to  quantify  spallation  resistance  of  soda-lime  glass 
and  identify  its  sensitivity  to  the  inclusion  of  non-linear  and  inelas¬ 
tic  effects. 

By  carrying  out  a  series  of  flyer  plate  impact  simulations  at 
increasing  impact  velocities,  the  critical  impact  velocity  beyond 
which  spallation  takes  place  has  been  determine  for  each  of  the 
three  aforementioned  material  model  cases.  Through  this  analysis 
procedure,  it  was  observed  that  the  critical  impact  velocity  is  ca. 
25  m/s  for  all  three  cases  considered.  This  finding  is  not  unex¬ 
pected  since:  (a)  due  to  a  low  value  of  the  fracture  strength  for 
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Fig.  10.  Temporal  evolution  of  the  target  back-face  particle  velocities  for  the  case  of 
flyer-plate  impact  velocities  of:  (a)  200  m/s;  and  (b)  900  m/s. 

soda-lime  glass,  one  should  expect  a  relatively  low  value  of  the 
critical  impact  velocity;  and  (b)  for  the  same  reason,  spallation  oc¬ 
curs  at  a  stress  level  where  material  nonlinearity  and  inelasticity 
effects  make  an  insignificant  contribution  to  the  overall  material 
response. 

The  finding  presented  above  shows  that  protection  performance 
of  soda-lime  glass  as  measured  by  the  critical  impact  velocity  is  not 
dependent  upon  material  nonlinearity  and  inelasticity.  However, 
one  may  explore  additional  aspects  of  the  impact  performance  of 
soda-lime  glass.  For  example,  momentum/kinetic  energy  carried 
by  the  propelled  principal  fragment  is  an  important  factor  relative 
to  the  threat  imposed  by  vehicle/structure-interior  fragment  re¬ 
lease.  In  other  words,  lower  values  of  the  principal  fragment  linear 
momentum  and  kinetic  energy  are  preferred  in  the  case  of  spalla¬ 
tion.  To  investigate  this  aspect  of  soda-lime  glass  impact  perfor¬ 
mance,  flyer-plate  experiments  are  simulated  at  an  impact 
velocity  of  900  m/s.  The  results  of  these  simulations  for  the  three 
material  model  cases  are  displayed  in  Fig.  lla-c,  respectively. 
Examination  of  the  results  displayed  in  these  figures  show  that 
selection  of  material  model  affects  both  the  size  and  the  velocity 
of  the  principal  fragment.  Specifically,  for  the  original,  nonlinear 
elastic,  and  nonlinear  elastic/inelastic  soda-lime  glass  material 
models  the  linear  momenta  (scaled  by  the  flyer-plate  initial 
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(a) 

Spall  Plane 

900m/s 

(b) 

Spall  Plane  _  ^ 

832m/s 

(c) 

Spall  Plane  _ _ _ 

825m/s 

Fig.  11.  Particle  velocity  field  plots  at  a  time  shortly  after  spallation  for  the 
following  material  model  cases  for  soda-linre  glass:  (a)  the  original  formulation 
[4,5];  (b)  the  nonlinear  elastic  formulation;  and  (c)  the  nonlinear  elastic/inelastic 
formulation. 


momentum)  are  approximately  1.0,  0.88,  and  0.90,  respectively. 
The  corresponding  scaled  kinetic  energy  of  the  principal  fragment 
are  approximately  1.0,  0.78,  and  0.77,  respectively.  These  findings 
suggest  that  nonlinear  elasticity  and  inelastic  deformation  provide 
beneficial  effects  in  reducing  the  linear  momentum  and  kinetic 
energy  of  the  principal  fragment.  It  should  be  noted,  that  the 
observed  reduction  in  the  principal  fragment  linear  momentum 
does  not  violate  the  principle  of  linear-momentum  conservation 
since  additional  momentum  remains  in  the  flyer-plate/target.  Also, 
at  least  a  portion  of  the  reduced  kinetic  energy  of  the  principal  frag¬ 
ment  is  associated  with  the  permanent  densification  of  a  portion  of 
the  soda-lime  glass  material  adjacent  to  the  flyer-plate  impact 
interface. 

The  findings  reported  above  can  be  directly  used  during  design 
of  mass-efficient  spallation-resistant  transparent-armor  systems. 
In  our  ongoing  work,  transparent  armor  design  guidelines  are 
being  developed  using  the  concept  of  material  selection  charts. 
The  foundation  for  this  approach  is  analogous  to  the  one  developed 
in  our  recent  work  [17], 

It  should  be  noted  that  experimental  investigation  of  fracture 
phenomena  (like  spall  fracture)  under  dynamic  loading  conditions 
in  brittle  materials  is  severely  limited  due  to  the  accompanying 
massive  destruction  of  the  test  samples.  Furthermore,  since  spall 
fracture  often  controls  ballistic  performance  of  transparent  armor 
systems,  it  is  usually  considered  to  be  a  sensitive  subject  matter. 
These  are  perhaps  the  two  main  reasons  for  paucity  of  relevant 
information  in  open  literature.  Typically,  reported  spall  fracture 
studies  involve  the  use  of  high-power  lasers  e.g.  [32]  to  generate 
shocks  (more  precisely,  short-duration  shock  pulses)  rather  than 
the  use  of  ballistic  impact  tests  e.g.  [33].  Since  a  most  comprehen¬ 
sive  set  of  the  open  literature  relevant  experimental  results  identi¬ 
fied  in  the  present  work  are  those  reported  in  Ref.  [32],  an  attempt 
is  made  in  the  remainder  of  this  manuscript  to  correlate  these 
experimental  results  with  the  present  computational  results.  This 
was  done  while  recognizing  that  the  results  in  Ref.  [32]  were  gen¬ 
erated  under  shock  pulse  loading  condition  while  the  computa¬ 
tional  analysis  carried  out  in  the  present  work  involved  sustained 
shock  loading.  Comparison  between  the  experimental  and  the 
computational  results  was  done  with  respect  to  the  effect  of:  (a) 
shock  strength;  and  (b)  target  plate  thickness  on  the  size  of  the 
damage  zone  and  the  final  size  of  the  principal  fragment. 

Experimental  results  reported  in  Ref.  [32]  showed  that  an 
increase  in  the  shock  strength  by  100%  led  to  a  60-70%  increase 
in  the  principal  spall  fragment  size.  The  present  computational 
results  show  that,  the  same  increase  in  shock  strength  results  in 
a  ca.  70%  increase  in  the  principal  spall  fragment  size.  As  far  as 
the  effect  of  target  plate  thickness  is  concerned,  experiments  in 
Ref.  [32]  show  that  an  increase  of  this  quantity  by  200%  resulted 


in  a  decrease  of  the  final  size  of  the  principal  fragment  by  ca. 
10%.  The  present  computational  results,  on  the  other  hand,  show 
that  the  principal  spall-fragment  size  is  independent  of  the  target 
thickness  but  rather  scales  directly  with  the  flyer  plate  (projectile) 
thickness.  It  should  be  noted  here  that  the  effect  of  the  target  thick¬ 
ness  observed  in  Ref.  [32]  is  mainly  associated  with  the  extent  of 
change  in  the  shock  pulse  profile  during  its  propagation  through 
the  target.  On  the  other  hand,  sustained  shocks  used  in  the  present 
computational  work  are  not  affected  by  the  target-plate  thickness. 
Based  on  the  aforementioned  outcome  of  the  experiment/compu¬ 
tation  comparison,  it  can  be  concluded  that  the  present  model 
yields  results  which  are  generally  consistent  with  the  open  litera¬ 
ture-reported  experimental  results. 

6.  Conclusions 

Based  on  the  material-model  development  procedure  utilized 
and  the  results  of  the  subsequent  computational  analyses,  the 
following  main  summary  remarks  and  conclusions  can  be  drawn: 

1.  The  effect  of  material  nonlinearity  and  inelastic  behavior  on  the 
dynamic  response  (including  spallation)  of  soda-lime  glass  is 
studied  under  symmetric  flyer-plate  loading  conditions  using 
computational  methods  and  tools. 

2.  The  nonlinear  elastic  and  inelastic  effects  are  incorporated  into 
the  high  strain-rate,  high-pressure,  large-strain  material  model 
for  soda-lime  ballistic  glass  recently  proposed  by  the  authors. 

3.  The  flyer-plate  impact  simulation  results  revealed  that  inclu¬ 
sion  of  nonlinear  elastic  and  inelastic  effects  do  not  measurably 
affect  spall  resistance  (as  measured  by  a  minimum  flyer-plate 
velocity  resulting  in  spallation). 

4.  However,  these  phenomena  are  found  to  yield  beneficial  effects 
associated  with  linear-momentum/kinetic  energy  reduction 
effects  in  the  spall  fragments. 
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